﻿/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms which obtain the physical parameters of Mars
Created: PJN / 04-01-2004
History: None

Copyright (c) 2003 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/

using System;

namespace AAPlus
{

    public class PhysicalMarsDetails
    {
        public double DE;
        public double DS;
        public double w;
        public double P;
        public double X;
        public double k;
        public double q;
        public double d;

    }

    public static class PhysicalMars
    {

        public static PhysicalMarsDetails Calculate(double JD)
        {
            //What will be the return value
            PhysicalMarsDetails details = new PhysicalMarsDetails();

            //Step 1
            double T = (JD - 2451545) / 36525;
            double Lambda0 = 352.9065 + 1.17330 * T;
            double Lambda0rad = CoordinateTransformation.DegreesToRadians(Lambda0);
            double Beta0 = 63.2818 - 0.00394 * T;
            double Beta0rad = CoordinateTransformation.DegreesToRadians(Beta0);

            //Step 2
            double l0 = Earth.EclipticLongitude(JD);
            double l0rad = CoordinateTransformation.DegreesToRadians(l0);
            double b0 = Earth.EclipticLatitude(JD);
            double b0rad = CoordinateTransformation.DegreesToRadians(b0);
            double R = Earth.RadiusVector(JD);

            double PreviousLightTravelTime = 0;
            double LightTravelTime = 0;
            double x = 0;
            double y = 0;
            double z = 0;
            bool bIterate = true;
            double DELTA = 0;
            double l = 0;
            double lrad = 0;
            double b = 0;
            double brad = 0;
            double r = 0;
            while (bIterate)
            {
                double JD2 = JD - LightTravelTime;

                //Step 3
                l = Mars.EclipticLongitude(JD2);
                lrad = CoordinateTransformation.DegreesToRadians(l);
                b = Mars.EclipticLatitude(JD2);
                brad = CoordinateTransformation.DegreesToRadians(b);
                r = Mars.RadiusVector(JD2);

                //Step 4
                x = r * Math.Cos(brad) * Math.Cos(lrad) - R * Math.Cos(l0rad);
                y = r * Math.Cos(brad) * Math.Sin(lrad) - R * Math.Sin(l0rad);
                z = r * Math.Sin(brad) - R * Math.Sin(b0rad);
                DELTA = Math.Sqrt(x * x + y * y + z * z);
                LightTravelTime = Elliptical.DistanceToLightTime(DELTA);

                //Prepare for the next loop around
                bIterate = (Math.Abs(LightTravelTime - PreviousLightTravelTime) > 2E-6); //2E-6 correponds to 0.17 of a second
                if (bIterate)
                    PreviousLightTravelTime = LightTravelTime;
            }

            //Step 5
            double lambdarad = Math.Atan2(y, x);
            double lambda = CoordinateTransformation.RadiansToDegrees(lambdarad);
            double betarad = Math.Atan2(z, Math.Sqrt(x * x + y * y));
            double beta = CoordinateTransformation.RadiansToDegrees(betarad);

            //Step 6
            details.DE = CoordinateTransformation.RadiansToDegrees(Math.Asin(-Math.Sin(Beta0rad) * Math.Sin(betarad) - Math.Cos(Beta0rad) * Math.Cos(betarad) * Math.Cos(Lambda0rad - lambdarad)));

            //Step 7
            double N = 49.5581 + 0.7721 * T;
            double Nrad = CoordinateTransformation.DegreesToRadians(N);

            double ldash = l - 0.00697 / r;
            double ldashrad = CoordinateTransformation.DegreesToRadians(ldash);
            double bdash = b - 0.000225 * (Math.Cos(lrad - Nrad) / r);
            double bdashrad = CoordinateTransformation.DegreesToRadians(bdash);

            //Step 8
            details.DS = CoordinateTransformation.RadiansToDegrees(Math.Asin(-Math.Sin(Beta0rad) * Math.Sin(bdashrad) - Math.Cos(Beta0rad) * Math.Cos(bdashrad) * Math.Cos(Lambda0rad - ldashrad)));

            //Step 9
            double W = CoordinateTransformation.MapTo0To360Range(11.504 + 350.89200025 * (JD - LightTravelTime - 2433282.5));

            //Step 10
            double e0 = Nutation.MeanObliquityOfEcliptic(JD);
            double e0rad = CoordinateTransformation.DegreesToRadians(e0);
            Coordinate2D PoleEquatorial = CoordinateTransformation.Ecliptic2Equatorial(Lambda0, Beta0, e0);
            double alpha0rad = CoordinateTransformation.HoursToRadians(PoleEquatorial.X);
            double delta0rad = CoordinateTransformation.DegreesToRadians(PoleEquatorial.Y);

            //Step 11
            double u = y * Math.Cos(e0rad) - z * Math.Sin(e0rad);
            double v = y * Math.Sin(e0rad) + z * Math.Cos(e0rad);
            double alpharad = Math.Atan2(u, x);
            double alpha = CoordinateTransformation.RadiansToHours(alpharad);
            double deltarad = Math.Atan2(v, Math.Sqrt(x * x + u * u));
            double delta = CoordinateTransformation.RadiansToDegrees(deltarad);
            double xi = Math.Atan2(Math.Sin(delta0rad) * Math.Cos(deltarad) * Math.Cos(alpha0rad - alpharad) - Math.Sin(deltarad) * Math.Cos(delta0rad), Math.Cos(deltarad) * Math.Sin(alpha0rad - alpharad));

            //Step 12
            details.w = CoordinateTransformation.MapTo0To360Range(W - CoordinateTransformation.RadiansToDegrees(xi));

            //Step 13
            double NutationInLongitude = Nutation.NutationInLongitude(JD);
            double NutationInObliquity = Nutation.NutationInObliquity(JD);

            //Step 14
            lambda += 0.005693 * Math.Cos(l0rad - lambdarad) / Math.Cos(betarad);
            beta += 0.005693 * Math.Sin(l0rad - lambdarad) * Math.Sin(betarad);

            //Step 15
            Lambda0 += NutationInLongitude / 3600;
            Lambda0rad = CoordinateTransformation.DegreesToRadians(Lambda0);
            lambda += NutationInLongitude / 3600;
            lambdarad = CoordinateTransformation.DegreesToRadians(lambda);
            e0 += NutationInObliquity / 3600;
            e0rad = CoordinateTransformation.DegreesToRadians(e0rad);

            //Step 16
            Coordinate2D ApparentPoleEquatorial = CoordinateTransformation.Ecliptic2Equatorial(Lambda0, Beta0, e0);
            double alpha0dash = CoordinateTransformation.HoursToRadians(ApparentPoleEquatorial.X);
            double delta0dash = CoordinateTransformation.DegreesToRadians(ApparentPoleEquatorial.Y);
            Coordinate2D ApparentMars = CoordinateTransformation.Ecliptic2Equatorial(lambda, beta, e0);
            double alphadash = CoordinateTransformation.HoursToRadians(ApparentMars.X);
            double deltadash = CoordinateTransformation.DegreesToRadians(ApparentMars.Y);

            //Step 17
            details.P = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(Math.Atan2(Math.Cos(delta0dash) * Math.Sin(alpha0dash - alphadash), Math.Sin(delta0dash) * Math.Cos(deltadash) - Math.Cos(delta0dash) * Math.Sin(deltadash) * Math.Cos(alpha0dash - alphadash))));

            //Step 18
            double SunLambda = Sun.GeometricEclipticLongitude(JD);
            double SunBeta = Sun.GeometricEclipticLatitude(JD);
            Coordinate2D SunEquatorial = CoordinateTransformation.Ecliptic2Equatorial(SunLambda, SunBeta, e0);
            details.X = MoonIlluminatedFraction.PositionAngle(SunEquatorial.X, SunEquatorial.Y, alpha, delta);

            //Step 19
            details.d = 9.36 / DELTA;
            details.k = IlluminatedFraction.ConvertToIlluminatedFraction(r, R, DELTA);
            details.q = (1 - details.k) * details.d;

            return details;
        }

    }
}
